R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(spatstat)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
install.packages("tidyverse")
library(tmap)
library(sf)
library(geojsonio)
library(tmaptools)
library(tidyverse)
library(rgdal)
library(raster)
library(SciViews)
library(base)
library(broom)
library(spdep)

####Basemap
library(tmap)
BoroughMap <- geojson_read("https://raw.githubusercontent.com/ft-interactive/geo-data/master/uk/london/london-wards-2014.geojson", what = "sp")
BNG = "+init=epsg:27700"
BoroughMapBNG <- spTransform(BoroughMap,BNG)
WGS = "+init=epsg:4326"
tmap_mode("view")
tm_shape(BoroughMapBNG) +tm_polygons(col = NA, alpha = 0.5)

#aggregate city of london 
library(sf)
library(dplyr)
library(rgdal)
BoroughMapSF <- st_as_sf(BoroughMapBNG)
library(rgeos)
CityofLondon<- gUnionCascaded(BoroughMapBNG[c(630:654), ])
plot(CityofLondon)
CityofLondonSF <- st_as_sf(CityofLondon)
CityofLondonSF$gss_code_ward<-"E05000001"  ##add same columns in city of london, for aggregating with other borough
CityofLondonSF$gss_code_borough <- "E09000001"
CityofLondonSF$borough <-"City Of London"
CityofLondonSF$ward <- "City Of London"
BoroughMapSFnew <- BoroughMapSF[-c(630:654), ]
LondonWardSF <- rbind(BoroughMapSFnew,CityofLondonSF) ##Attention New not new
LondonWard <- as(LondonWardSF,"Spatial")
rownames(LondonWardSF) <- 1:nrow(LondonWardSF)
tmap_mode("view")
tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)#+tm_text("ward",size=0.3, col="black") 

####population density (using https://data.london.gov.uk/dataset/land-area-and-population-density-ward-and-borough 2016)
library(tidyverse)
Density <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/population_density_2016_final.csv",na = "n/a")
#normalization population_per_square_km
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
density.normalization <- normalization(Density$`population density`)
Density$Density.normalization <- density.normalization

#join the data to basemap P3P27
LondonWard.join.popdens <- LondonWardSF%>% left_join(Density,by=c("gss_code_ward"="New Code"))
tmap_mode("view")
qtm(LondonWard.join.popdens,fill="Density.normalization")

### Public transport accessibility level(PTAL)
PTAL <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Ward2014 Avpublic transport accessibility level2015.csv",na = "n/a")
#normalization PTAL
normalization<-function(x){
return((x-min(x))/(max(x)-min(x)))}
PTAL.normalization <- normalization(PTAL$AvPTAI2015)
PTAL$PTAL.normalization <- PTAL.normalization
#join the data to basemap P3P27
LondonWard.join.PTAL <- LondonWardSF%>% left_join(PTAL,by=c("gss_code_ward"="Ward Code"))
qtm(LondonWard.join.PTAL,fill="PTAL.normalization")
qtm(LondonWard.join.PTAL,fill="AvPTAI2015")

###Public open space accessibility level(POSAL)
POSAL <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/access_public_open_space_ward.csv",na = "n/a")
POSAL$POSAL.normalization <- normalization(POSAL$`Open Space`)
LondonWard.join.POSAL <- LondonWardSF%>% left_join(POSAL,by=c("gss_code_ward"="WD13CD"))
#qtm(LondonWard.join.POSAL,fill="POSAL.normalization")
#combine two accessibility sub-index
Access.normalization<- merge(POSAL,PTAL, by.x="WD13CD",by.y="Ward Code",all=TRUE)
#remove useless column and rows
Access.normalization[,c(4:8,10:11)] <- NULL
#calculate the mean of two accessibility level
Access.normalization$access.mean <- rowMeans(Access.normalization[c('PTAL.normalization', 'POSAL.normalization')], na.rm=TRUE)

node

### Road intersections
node <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/Road Networkoproad_essh_gb/new_road_node.shp",layer="new_road_node")
BNG = "+init=epsg:27700"
nodeBNG <- spTransform(node, BNG)
summary(node)
##qtm(nodeBNG)
nodeBNGSF <- st_as_sf(nodeBNG)
#select the node type"junction"
nodeBNGSF <- nodeBNGSF[which(nodeBNGSF$formOfNode=="junction"),]
nodeBNGSP <- as(nodeBNGSF, "Spatial")
head(node)
#select node point in borough
nodeBNGSP <- remove.duplicates(nodeBNGSP)
nodeBNGSP.final<- nodeBNGSP[LondonWard,]
# tmap view
#tmap_mode("view")
#tm_shape(LondonWard) +
#  tm_polygons(col = NA, alpha = 0.5) +
#tm_shape(nodeBNGSP.final) +
#  tm_dots(col = "blue")

# count number of point in each ward 
library(GISTools)
poly.counts(nodeBNGSP.final, LondonWard) -> nodecount
##setNames(nodecount, LondonWard@data$nodecount)
London.node <- LondonWardSF
London.node$nodecount <- nodecount   #add number of point to each ward (London.node$nodecount)
#combine the ward area(sq km)
London.node$Square_Kilometres <- LondonWard.join.popdens$`Square Kilometres`
#calculate point/area--point per square kilometre
London.node$PointPerSK <-London.node$nodecount/London.node$Square_Kilometres
qtm(London.node,fill="PointPerSK")
#join with LondonWardSF
London.nodeSF <- append_data(LondonWardSF,London.node, key.shp = "code", key.data = "gss_code_ward", ignore.duplicates = TRUE)

economic vitality

###economic vitality
#CLASS 1  house price
#HousePrice <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/house_price_mean2016.csv",na = "n/a")
#HousePrice.normalization <- normalization(HousePrice$Value)
#HousePrice$HousePrice.normalization <- HousePrice.normalization
#CLASS 2 economic activity
EcoAct <- read_csv("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/economic_activity_ward.csv",na = "n/a")
EcoAct.normalization <- normalization(EcoAct$`Economically active: Total`)
EcoAct$EcoAct.normalization <- EcoAct.normalization

land use mix

### Land use mix
#https://en.wikipedia.org/wiki/Planning_use_classes_in_England 
# six classes--
landcover <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/gis_osm_landuse_a_free_1.shp",layer="gis_osm_landuse_a_free_1")
landcoverBNG <- spTransform(landcover, BNG)
landcoverBNGSF <- st_as_sf(landcoverBNG)

AOIS <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/gis_osm_pois_a_free_1.shp",layer="gis_osm_pois_a_free_1")
AOISBNG <- spTransform(AOIS, BNG)
AOISBNGSF <- st_as_sf(AOISBNG)

#find missing data in each column in landcoverBNGSF
colSums(is.na(landcoverBNGSF))
#get the classes of landcover
summary(landcoverBNGSF)
#CLASS1 greenspace
greenspaceBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="park" |  landcoverBNGSF$fclass =="forest"|landcoverBNGSF$fclass =="nature_reserve"|landcoverBNGSF$fclass =="grass"),] 
greenspaceBNGSP <- as(greenspaceBNGSF,"Spatial")
#CLASS2 industrial and commercial
industrialBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="industrial"| landcoverBNGSF$fclass =="commercial"),]
industrialBNGSP <- as(industrialBNGSF,"Spatial")
#CLASS 3 residential
residentialBNGSF <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="residential" ),]
residentialBNGSP <- as(residentialBNGSF,"Spatial")
#CLASS 4 recreation-leisure
recreationBNGSF1 <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="recreation_ground"| landcoverBNGSF$fclass =="allotments"),]
recreationBNGSF2 <- AOISBNGSF[ which(AOISBNGSF$fclass =="attraction"|AOISBNGSF$fclass =="theatre"|AOISBNGSF$fclass =="arts_centre"),]
recreationBNGSF <- rbind(recreationBNGSF2,recreationBNGSF1)
recreationBNGSP <- as(recreationBNGSF,"Spatial")


#CLASS 5 retail
retailBNGSF1 <- landcoverBNGSF[ which( landcoverBNGSF$fclass =="retail" ),]
retailBNGSF2 <- AOISBNGSF[ which( AOISBNGSF$fclass =="retail" ),]
retialBNGSF <- rbind(retailBNGSF2,retailBNGSF1)
retailBNGSP <- as(retailBNGSF,"Spatial")

#CLASS 6--community services
communityBNGSF <- AOISBNGSF[ which(AOISBNGSF$fclass =="school"|AOISBNGSF$fclass =="hospital"|AOISBNGSF$fclass =="nursing_home"|AOISBNGSF$fclass =="university"  ),]
communityBNGSP <- as(communityBNGSF,"Spatial")

#CLASS 7--transport
transport <- readOGR("E:/UCL/005-GI System&Science/GIS-Assessment 3/R assessment3/greater-london-latest-free.shp/road_polygon_merge.shp",layer="road_polygon_merge")
transportBNG <- spTransform(transport, BNG)
transportBNG <- transportBNG[transportBNG@data$AREA >10000,] 

# Assign the land use area for each subset. Here we first use the most widespread residential as example 
library(raster)
r <- raster(ncols=400, nrows=400) #generate raster size(quantity)
extent(r) <- extent(residentialBNGSP)  #Important! Assign the extent of raster to cover the same extents of the polygon
residential.Raster <- rasterize(residentialBNGSP,r,background=NA) #convert polygon to raster
#tmap_mode("view")
#qtm(residential.Raster)+tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.5)
#extract rasters from polygons#(extract() cannot use method=bilinear in raster-in-polygon)
residential.extract <- raster::extract(residential.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
#remove rows with NA
row.has.na <- apply(residential.extract, 1, function(x){any(is.na(x))})
residential.extract2 <-residential.extract[!row.has.na,]
#merge the pixel in same ward(aggregate the frequency of the same data in one)
library(plyr)
residential.extract.sum <- factor(residential.extract2$ID)  
residential.extract.sum <- table(residential.extract2$ID)
residential.extract.wardsum <- as.data.frame(residential.extract.sum)
#create new order number for LondonWardSF, for merge the residential.extract and LondonWardSF by order
LondonWard2SF <- LondonWardSF
LondonWard2SF$Var1 <- 1:630
LondonWard2SF[,'Var1']<-factor(LondonWard2SF[,'Var1'])
#merge two dataframe
residential.extract.final <- merge(residential.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
#residential.extract.final [!duplicated(residential.extract.final ), ]
#replace NA values with 0
residential.extract.final[is.na(residential.extract.final)] <- 0
#rename the column
colnames(residential.extract.final)[2] <- "Freq.residential"
#raster cell size
#we have to convert raster to polygon because of raster(). This function only compute the pixel area in longitude/latitude coordiante system
residential.repolygon <- rasterToPolygons(residential.Raster, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=TRUE)
#qtm(residential.repolygon)
residential.repolygonSF <- st_as_sf(residential.repolygon) 
residential.repolygonSF$area <- st_area(residential.repolygonSF)
#so the single pixel area is 28999.54 [m^2](0.029 km^2)

land use mix functions

##based on the details above, we create functions for land use analysis
function.raster <- function(landsp){
  extent(r) <- extent(landsp)
  landsp.Raster <- rasterize(landsp,r,background=NA)
  return(landsp.Raster)
}

function.extract<- function(landsp.Raster){
  landsp.extract <- raster::extract(landsp.Raster,LondonWard,df=TRUE, weights =FALSE, na.rm = TRUE)
  row.has.na <- apply(landsp.extract, 1, function(x){any(is.na(x))})
  landsp.extract2 <-landsp.extract[!row.has.na,]
  landsp.extract.sum <- factor(landsp.extract2$ID)  
  landsp.extract.sum <- table(landsp.extract2$ID)
  landsp.extract.wardsum <- as.data.frame(landsp.extract.sum)
  landsp.extract.final <- merge(landsp.extract.wardsum, LondonWard2SF, by="Var1",all=TRUE)
  landsp.extract.final[is.na(landsp.extract.final)] <- 0
  return(landsp.extract.final)
}

#green space raster
greenspace.Raster <- function.raster(greenspaceBNGSP)
greenspace.extract.final <- funtion.extract(greenspace.Raster)
colnames(greenspace.extract.final)[2] <- "Freq.greenspace"

#industrial and commercial
industrial.Raster <- function.raster(industrialBNGSP)
industrial.extract.final <- function.extract(industrial.Raster)
colnames(industrial.extract.final)[2] <- "Freq.industrial"
#recreation
recreation.Raster <- function.raster(recreationBNGSP)
recreation.extract.final <- function.extract(recreation.Raster)
colnames(recreation.extract.final)[2] <- "Freq.recreation"
#retail
retail.Raster <- function.raster(recreationBNGSP)
retail.extract.final <- function.extract(retail.Raster)
colnames(retail.extract.final)[2] <- "Freq.retail"
#community
community.Raster <- function.raster(communityBNGSP)
community.extract.final <- function.extract(community.Raster)
colnames(community.extract.final)[2] <- "Freq.community"
#transport
transport.Raster <- function.raster(transportBNG)
transport.extract.final <- function.extract(transport.Raster)
colnames(transport.extract.final)[2] <- "Freq.transport"

tmap_mode("view")
tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.2)+qtm(recreation.Raster)
#tm_shape(LondonWard) +tm_polygons(col = NA, alpha = 0.2)+qtm(residential.Raster)

#extract useful columns from multiple dataframe
landuse.dataframe <- residential.extract.final
landuse.dataframe <- merge(landuse.dataframe,greenspace.extract.final, by="Var1",all=TRUE)
landuse.dataframe <- merge(landuse.dataframe,industrial.extract.final, by="Var1",all=TRUE)
landuse.dataframe <- merge(landuse.dataframe,recreation.extract.final, by="Var1",all=TRUE)
landuse.dataframe <- merge(landuse.dataframe,retail.extract.final, by="Var1",all=TRUE)
landuse.dataframe <- merge(landuse.dataframe,community.extract.final, by="Var1",all=TRUE)
landuse.dataframe <- merge(landuse.dataframe,transport.extract.final, by="Var1",all=TRUE)
landuse.dataframe[,c(9:13,15:19,21:25,27:31,33:37,39:43)] <- NULL #remove useless column
colnames(landuse.dataframe)[2] <- "residential.freq" #rename column
head(landuse.dataframe)

index.dataframe[is.na(index.dataframe)] <- 0 
###compute land use mix index
index.dataframe$MixIndex <- (abs(rowSums(index.dataframe[,8:14])))/ln(7)
#MixIndex normalisation
index.dataframe$MixIndex.normalisation <- normalization(MixIndex )

Error in normalization(MixIndex) : object 'MixIndex' not found

qtm(residential.Raster)
library(tmap)
tmap_mode("view")
tm_shape(LondonWard) +
  tm_polygons(col = NA, alpha = 0.5) +
tm_shape(residentialBNGSP) +
  tm_polygons(col = "blue")

Test for the relationship

Linear regression model

###Since we have done the most important part of data manipulation, next we will explore the relationship between intersection and vitality
#histogram of two variables
qplot(PointPerSK, data = London.node, geom = "histogram")

#ggplot(London.node, aes(x="PointPerSK", y=PointPerSK))+geom_boxplot()#(this is another code for boxplot in ggplot2)
#qplot(y=London.node$PointPerSK, x= 0.1, geom = "boxplot")  
qplot(Index, data = index.dataframe, geom = "histogram")

#qplot(y= index.dataframe$Index, x= 0.1, geom = "boxplot")  
#the variables look normally distributed,but we should consider the outliers
plot(London.node$PointPerSK,index.dataframeSF$Index,col="blue")

#after scatter plot, we should remove "Lady Margaret" as its unusual road intersection density(maybe due to the residential road intersections are included)
node.index.frame <-merge(index.dataframe,London.node, by.x="gss_code_ward",by.y="gss_code_ward",all=TRUE) 
#remove outlier--Lady Margaret(E05000181)
#node.index.frame <- node.index.frame[node.index.frame$PointPerSK<300,]
#linear regression plot
qplot(PointPerSK,Index,data=node.index.frame,geom = "point") + stat_smooth(method="lm", se=FALSE, size=1)

#fit the linear model
model <- lm(PointPerSK ~ Index, data =node.index.frame)
#residual dataframe
library(broom)
model_res <- tidy(model)
summary(model)

Call:
lm(formula = PointPerSK ~ Index, data = node.index.frame)

Residuals:
    Min      1Q  Median      3Q     Max 
-80.543 -20.260  -5.167  13.257 198.842 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -12.822      6.281  -2.041   0.0416 *  
Index        246.958     14.180  17.416   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.48 on 628 degrees of freedom
Multiple R-squared:  0.3257,    Adjusted R-squared:  0.3246 
F-statistic: 303.3 on 1 and 628 DF,  p-value: < 2.2e-16

plot(model)

Fitting a Quadratic Model

node.index.frame$PointPerSK2 <- (node.index.frame$PointPerSK)^2
quadratic.model <-lm(Index ~ PointPerSK + (PointPerSK2),data =node.index.frame)
plot(quadratic.model )

summary(quadratic.model)

Call:
lm(formula = Index ~ PointPerSK + (PointPerSK2), data = node.index.frame)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.186415 -0.045824  0.002214  0.043010  0.277679 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.861e-01  1.235e-02   15.07   <2e-16 ***
PointPerSK   3.927e-03  2.265e-04   17.34   <2e-16 ***
PointPerSK2 -1.175e-05  9.728e-07  -12.08   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06766 on 627 degrees of freedom
Multiple R-squared:  0.4529,    Adjusted R-squared:  0.4512 
F-statistic: 259.6 on 2 and 627 DF,  p-value: < 2.2e-16

#node.index.frame$Index<- 1.861e-01+(3.927e-03)*(node.index.frame$PointPerSK)+(-1.175e-05)*(node.index.frame$pointPerSK2)
#plot(node.index.frame$PointPerSK,node.index.frame$Index,type='l',col='navy',main='Nonlinear relationship',lwd=3)
PointPerSKvalues <- seq(0,300, 0.01)
predictedcounts <- predict(quadratic.model,list(PointPerSK=PointPerSKvalues, PointPerSK2=PointPerSKvalues^2))
plot(node.index.frame$PointPerSK, node.index.frame$Index , pch=16, xlab = "PointPerSK", ylab = "Index", cex.lab = 0.8, col = "black")
lines(PointPerSKvalues, predictedcounts, col = "blue", lwd =2)

The coefficient (estimate) in the summary() table shows that(on average) for a 1 unit (day) increase in unauthorised absence from school, there is a reduction in the average GCSE point score of -29.874.

The p-values for the intercept and the coefficient are highly statistically significant (<0.001) so we can rely on the relationship that is being observed.

The adjusted R-squared statistic is 0.38, which tells us that 38% of the variation in GCSE scores across Wards in London can be explained by variation in unauthorised absence from school (which is quite a lot for a single variable).

Interrogating the last graph in plot(model1)which is a scatter plot of fitted values (the model estimates achieved by plugging the values and coefficients back into the regression equation) against standardised residuals, we can see no apparent patterns in the cloud of points, which suggests the model has not violated any important assumptions.

Test for spatial patterns

we shoud test that is there any spatial clustering of residuals.If these places cluster in space, then there might be some unobserved underlying factor causing this. This is important is it means that the assumption of independence that regression models rely upon might be violated. Now in our case here, there is not a clear violation of spatial independence, but it is certainly hinted at. test for spatial patterns (spatial autocorrelation) using the Moran’s I statistic

### Test for spatial patterns (spatial autocorrelation) using the Moran's I statistic
#copy a new LondonWardSF
LondonWardSF.pattern <- LondonWardSF
LondonWardSF.pattern$model_resids <- model$residuals
LondonWardSF.pattern <- LondonWardSF.pattern[!(is.na(LondonWardSF.pattern$model_resids) ),]
LondonWardSP.pattern <- as(LondonWardSF.pattern,"Spatial")
tmap_mode("view")

tmap mode set to interactive viewing

qtm(LondonWardSF.pattern,fill='model_resids')

Variable "model_resids" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

#library(spdep)
#Calculate the centroids of London wards
centWard <- coordinates(LondonWardSP.pattern)
plot(centWard)

#Generate the spatial weights matrix,using binary matrix of queen's case neighbours
#create neighbours list
library(spdep)
neighbour <- poly2nb(LondonWardSP.pattern,queen=T)
#plot
plot(neighbour,coordinates(centWard),col="red")
plot(LondonWardSP.pattern,add=T)

#create spatial weights object from these weights
SpatialWeight <- nb2listw(neighbour,style="C",zero.policy=T)
#Run moran I test
moran.test(LondonWardSP.pattern@data$model_resids, SpatialWeight,zero.policy=T)

    Moran I test under randomisation

data:  LondonWardSP.pattern@data$model_resids  
weights: SpatialWeight  n reduced by no-neighbour observations
  

Moran I statistic standard deviate = 9.827, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.2252039018     -0.0015923567      0.0005326405 

so here we have a statistically significant but relatively weak indication that there is some spatial clustering of residual values. A value of 0.07 (1 being perfect spatial autocorrelation, 0 being none at all) shows that there is no evidence that high residual values cluster near high values and low residual values cluster near lower values.

multi linear regression model

###For a deeper research, we will build a multiple linear regression model to find whether the urban vitality index has a relationship with urban vitality sub-factors. If so, which one is the strongest or weakest.
#
library(corrplot)
cormat <- cor(node.index.frame[,c(5:8,16)], use="complete.obs", method="pearson")
corrplot(cormat)

#rRun multiple regression model
model_multi <- lm(PointPerSK ~ MixIndex.normalization+Density.normalization + access.mean + EcoAct.normalization, data = node.index.frame)
summary(model_multi)

Call:
lm(formula = PointPerSK ~ MixIndex.normalization + Density.normalization + 
    access.mean + EcoAct.normalization, data = node.index.frame)

Residuals:
    Min      1Q  Median      3Q     Max 
-62.704 -15.816  -3.713  12.011 218.380 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)              26.098      5.856   4.457 9.86e-06 ***
MixIndex.normalization   16.078      7.086   2.269 0.023615 *  
Density.normalization   150.627      6.213  24.245  < 2e-16 ***
access.mean              26.792      7.021   3.816 0.000149 ***
EcoAct.normalization      8.127      6.392   1.271 0.204039    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.27 on 625 degrees of freedom
Multiple R-squared:  0.5611,    Adjusted R-squared:  0.5583 
F-statistic: 199.8 on 4 and 625 DF,  p-value: < 2.2e-16

summary(model_multi)

Call:
lm(formula = PointPerSK ~ MixIndex.normalization + Density.normalization + 
    access.mean + EcoAct.normalization, data = node.index.frame)

Residuals:
    Min      1Q  Median      3Q     Max 
-62.704 -15.816  -3.713  12.011 218.380 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)              26.098      5.856   4.457 9.86e-06 ***
MixIndex.normalization   16.078      7.086   2.269 0.023615 *  
Density.normalization   150.627      6.213  24.245  < 2e-16 ***
access.mean              26.792      7.021   3.816 0.000149 ***
EcoAct.normalization      8.127      6.392   1.271 0.204039    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.27 on 625 degrees of freedom
Multiple R-squared:  0.5611,    Adjusted R-squared:  0.5583 
F-statistic: 199.8 on 4 and 625 DF,  p-value: < 2.2e-16